function [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol)

  vP{1}=vectorP(maxBOs(1), aMinBE(1), aMaxBE(1));
  vP{2}=vectorP(maxBOs(2), aMinBE(2), aMaxBE(2));

  %Neccesary for speed-up
  vLambda=sort(vLambda);
  
  mPi=findPi2(CW,maxBOs,vM,N,vLambda,vP,step,errTol);

  %Indexes in the Results Vector
  indRV.lambda=  1;
  indRV.error=   2;
  indRV.S=       3;
  indRV.PDR=     4;
  indRV.pSend=   5;
  indRV.latency= 6;
  indRV.Pi=      7;
  indRV.Pt=      8;
  indRV.WRx=     9;
  indRV.WTx=    10;
  indRV.WIdle=  11;
  indRV.pReject=12;
  
  %Matrix with the estimated Results
  resCMatrix={};
  for i=1:length(vLambda)
    %Solve the model
    vPi=mPi(i,:);
    %Compute solution with a given mPi
    [vAuxPi, vS, vPt, vPDR, pSend, vLatency, vWRx, vWTx, vWIdle, pReject]=solveModel(CW, maxBOs, vM, N, vLambda(i), vP, vPi);

    resCVector{indRV.error}=   max(abs(vPi-vAuxPi));
    resCVector{indRV.S}=       vS;
    resCVector{indRV.PDR}=     vPDR;
    resCVector{indRV.pSend}=   pSend;
    resCVector{indRV.latency}= vLatency;
    resCVector{indRV.Pi}=      vPi;
    resCVector{indRV.Pt}=      vPt;
    resCVector{indRV.lambda}=  vLambda(i);
    resCVector{indRV.WRx}=     vWRx;
    resCVector{indRV.WTx}=     vWTx;
    resCVector{indRV.WIdle}=   vWIdle;
    resCVector{indRV.pReject}= pReject;

    resCMatrix=[resCMatrix; resCVector];
  end %for i=1:length(

end %function [

function mPi=findPi2(CW,maxBOs,vM,N,vLambda,vP,step,errTol)

  %Matrix with the estimated Results of pi
  mPi=[];
  tmpPi=1;

  % Sweep inter arrival rate in packets per packet duration 0 <= lambda <= 1
  for lambda=vLambda
    error=inf;
    %Empty the Solution Vector
    mPi=[mPi;0 0];
    % Guessing Pi. Probability that the channel is iddle
    for pi=tmpPi:-step:step
      %Prob that the channel is idle 'ind' consecutive slots
      %vPi(2) -> This assumption is true if and only if N1=N2=...=Nn, Eqs(5-6)
      vPi =[pi, pi*(N*pi-1+pi)/(N*pi)];
      %Solve the model
      [vAuxPi]=solveModel(CW, maxBOs, vM, N, lambda, vP, vPi);
      %Compare guessed values pi&pi_ii with new ones computed after solving the models
      tmpErr=abs(pi-vAuxPi(1));
      if(error>tmpErr)
        error=tmpErr;
        mPi(end,:)=vPi;
      end %if(error>
      if(error<errTol);  break; end
    end %for pi=vPi
    tmpPi=mPi(end,1);
    disp(['Lambda=', num2str(lambda), ', pi=', num2str(tmpPi), ', Error=', num2str(error)]);
  end %for lamda=

end


function [vAuxPi, vS, vPt, vPDR, pSend, vLatency, vWRx, vWTx, vWIdle, pReject]=solveModel(CW, maxBOs, vM, N, lambda, vP, vPi)
  %Maximum CW
  eCW=CW(end);

  %Default value in case of error
  vAuxPi=ones(1,eCW)*NaN;
  
  %Probability of transmission start Eq. (1)
  p=lambda/N;
  
  %Cond prob. that the channel is idle, knowing that it was idle in the previous 'ind-1' backoff slots. Eqs(5-6)
  vPi_=vPi(1);
  %This assumption is true if and only if N1=N2=...=Nn
  if(eCW>1)
    vPi_(2)=(N*vPi(1)-1+vPi(1))/(N*vPi(1));
  end
  %for i=3:eCW
  %  vPi_(i)=vPi(i)/vPi(i-1);
  %end

  %% Solve states for CAP for each type of nodes grouped by CW
  for i=1:length(CW)
    [IDLE, TX, BO, CS, tmpErr]=solveNode(p, vPi_,vP{i},maxBOs(i),CW(i));
    if(tmpErr),  return;  end
    %Summation of all possible channel states Eq(4). See Ramachandran Eq(11) for simplification details
    sumTmp=(1+TX*(N-1));
    %Probability that a Node 'ind' transmit in a generic backoff slot Eq(3)
    vPt(i)=vPi_(CW(i))*sum(CS(end,:))/sumTmp;
    if(nargout>=1)
      %Steady state probability to be in the TX state (Fraction of time spent in TX)
      vPtx(i)=N*TX/sumTmp;    
      %Steady state probability to be in the IDLE state (Fraction of time spent in IDLE)
      vPidle(i)=IDLE/sumTmp;    
      %Steady state probability to be in the BackOff, BOx, state (Fraction of time spent in BOx)
      vPbo(i)=sum(BO)/sumTmp;
      %Steady state probability to be in the Clear Channel Assessment (CCA), CSxy, state  (Fraction of time spent in CSxy)
      vPcs(i)=sum(sum(CS))/sumTmp;
      %Steady state probability to be in the first Clear Channel Assessment State (CCA), CSx1, state  (Fraction of time spent in CS1x)
      vPcs_1(i)=sum(CS(1,:))/sumTmp;
      
      %Prob that a node tryes to send a packet = Backoff Stages sending / Backoff stages the node has to be sending
      pSend(i)=vPtx(i)/(lambda*vPidle(i));
    end
    %Probability that a Node begin trans. given that it has sensed the channel idle 'ind' consecutive backoff slots Eq(7)
    vPt_(i)=vPt(i)./vPi(CW(i));
  end
    
  if(CW==[1 2])
    %Two types of nodes with CW=1,2

    %Alpha, no one (that is allowed to) begin a transmission Eq(8)
    vA(1)=(1-vPt_(1))^vM(1);
    vA(2)=vA(1)*(1-vPt_(2))^vM(2);
    
    %Beta, just one do it (mB(destination,source)) Eq(9)
    mB(1,1)=vM(1)*(1-vPt_(1))^(vM(1)-1)*vPt_(1);                      %From IDLE,BUSY to Success1
    mB(1,2)=mB(1,1)*(1-vPt_(2))^(vM(2));                              %From IDLE,IDLE to Success1

    mB(2,1)=0;                                                        %From IDLE,BUSY to Success2
    mB(2,2)=vM(2)*(1-vPt_(2))^(vM(2)-1)*vPt_(2)*(1-vPt_(1))^(vM(1));  %From IDLE,IDLE to Success2
  elseif(CW==[2 2])
    %Two types of nodes both with CW=2

    %Alpha, no one (that is allowed to) begin a transmission Eq(8)
    vA(1)=1;
    vA(2)=(1-vPt_(1))^vM(1)*(1-vPt_(2))^vM(2);

    %Beta, just one do it (mB(destination,source)) Eq(9)
    mB(1,1)=0;                                                        %From IDLE,BUSY to Success1
    mB(1,2)=vM(1)*(1-vPt_(1))^(vM(1)-1)*vPt_(1)*(1-vPt_(2))^(vM(2));  %From IDLE,IDLE to Success1

    mB(2,1)=0;                                                        %From IDLE,BUSY to Success2
    mB(2,2)=vM(2)*(1-vPt_(2))^(vM(2)-1)*vPt_(2)*(1-vPt_(1))^(vM(1));  %From IDLE,IDLE to Success2
  else
    error('Contention Window parameter nos supported, this solution is valid for CW=[1 2] or CW=[2 2]');
  end
  
  %Solve states for channel Markov Chain
  [IBUSY, SUCCESS, FAIL, tmpErr]=solveChannel(vA,mB,eCW);
  if(tmpErr),  return;  end 

  %% Compute vAuxPi from channel
  %Sum of dwell time in all channel states Eq(13)
  sumChSt =sum(IBUSY)+N*(sum(SUCCESS)+FAIL);

  %Prob that the channel is idle CW slots Eq(12)
  vPi_b=IBUSY(1:eCW-1)/sumChSt;
  
  %Cond prob. that the channel is idle, knowing that it was idle in the previous 'ind-1' backoff slots Eq(15)
  vAuxPi(eCW)=IBUSY(end)/sumChSt;

  for i=1:eCW-1
    %Prob that the channel is idle 'ind' consecutive slots
    vAuxPi(i)=vAuxPi(eCW)+sum(vPi_b(i:eCW-1));
  end
  
  %% Just vAuxPi??
  if(nargout==1)
    return;
  end
  
  %Throughput [Stot, S1, S2, ..., Sn] per class. Effective speed (KBps) = Sx/Mx* 250e3Kbps/(8b/B) Eq(16)
  vS=N*SUCCESS/sumChSt;
  vS=[sum(vS), vS];
  
  % Packet Delivery Ratio Type 2: #successfully received packets/#sent packets = sum(S)/sum(Mi*Ptxi) Eq(26)
  vPDR=[vS(1)/sum(vPtx.*vM), vS(2:end)./(vPtx.*vM)];
  
  %Latency in backoff slots Eq(27)
  vLatency=N*vM./vS(2:end).*(1-vPidle);
  
  %Power consumption Chipcon cc2420 802.15.4 compliant RF transceiver
  Nbeacon=2;              % Beacon length (in backoff slots)
  Nir=0.6;                % This is the duration of the switching period idle->receive (in backoff slots)
  BI=3072;                % Beacon Interval (in backoff slots)
  Widle=712e-3;           % Power consumed in idle state (mW)
  Wrx=35.28;              % Power consumed in RX state (mW)
  Wtx=31.32;              % Power consumed in TX state (mW)

  Pbeacon=Nbeacon/BI;     %Probability to receive a beacon -> Beacon Length/Beacon Interval length
  Pir=Nir*(1/BI+vPcs_1);  %Prob to be activating the radio Time to switch idle-receive/BI due to be in CSx1 or reciving Beacon
  
  %Power cosumed during reception Eq(20)
  vWRx=(vPcs+Pbeacon+Pir)*Wrx;
  %Power cosumed during transmission Eq(19)
  vWTx=vPtx*Wtx;
  %Power cosumed during idle period Eq(18)
  vWIdle=(vPidle-Pbeacon-Pir+vPbo)*Widle; 

  %Probability of throw away a new packet Eq(23)
  pReject=(1-vPidle);

end

